!-------------------------------------------------------------------------------------
! subroutine UpWind2ndOrder
!-------------------------------------------------------------------------------------

subroutine UpWind2ndOrder(  u,  dt, dx, a, j, uLeft,  uRight )

  implicit none

  integer j, i
  real u(j),  dt, dx, a, uLeft,  uRight
  real sigma 

  ! Vol. 2 page 502 of Numerical Computation of Internal and External Flows by Hirsch
  ! du/dt = -a du/dx
  ! For a >0

  ! sigma = a * dt / dx ,  sigma < 2 for stability

  ! u(i,n+1) = u(i) - (sigma/2)*[ 3*u(i) - 4*u(i-1) + u(i-2) ] + (sigma*sigma / 2 )*[ u(i) - 2*u(i-1) + u(i-2) ]

  !            all u's on RHS are at time step "n"

  if (j .eq. 1) return

  sigma = a * dt / dx

  !write(*,*) 'sigma=', sigma,'   j=',j  ; pause 'pause 01 in subroutine'
  !write(*,*)'u=', u     ; pause 'pause 02 in subroutine'

  if( sigma > 0.0 ) then  
     do i=j,3,-1
        u(i) = u(i) - (sigma/2)*( 3*u(i) - 4*u(i-1) + u(i-2) ) + &
             (sigma*sigma / 2 )*( u(i) - 2*u(i-1) + u(i-2) )
     enddo
     u(2) = u(2) - sigma*( u(2) - u(1) )!  first order for u(2)
     u(1) = uLeft
  else  
     do i=1, j-2
        u(i) = u(i) + (sigma/2)*( 3*u(i) - 4*u(i+1) + u(i+2) ) + &
             (sigma*sigma / 2 )*( u(i) - 2*u(i+1) + u(i+2) )
     enddo
     u(j-1) = u(j-1) - sigma*( u(j) - u(j-1) )!  first order for u(2)
     u(j) = uRight
  endif

  return

end subroutine UpWind2ndOrder


!-------------------------------------------------------------------------------------
! subroutine UpWind2ndOrder_Variable_a
!-------------------------------------------------------------------------------------

subroutine UpWind2ndOrder_Variable_a( u, dt, dx, aVar, j, uLeft,  uRight, iPeriodic )

  ! as per communication from B. Wendroff on Jan 10, 2011, for:  du/dt =  - (d/dx)[ aVar(x) u ]
  ! see Word doc:   Advection in non uniform flow

  implicit none

  integer j, i, iPeriodic
  real u(j), dt, dx, aVar(j), uLeft,  uRight
  real sigma
  real    u0,    uM1,    uJp1,    uJp2
  real aVar0, aVarM1, aVarJp1, aVarJp2

  ! Vol. 2 page 502 of Numerical Computation of Internal and External Flows by Hirsch
  ! du/dt = -aVar du/dx for aVAr constant
  ! sigma = aVar_max * dt / dx ,  sigma < 2 for stability if aVar is uniform

  ! u(i,n+1) = u(i) - (sigma/2)*[ 3*u(i) - 4*u(i-1) + u(i-2) ] + (sigma*sigma / 2 )*[ u(i) - 2*u(i-1) + u(i-2) ]

  !            all u's on RHS are at time step "n"

  if (j .le. 2) return

  sigma =     dt / dx

  ! needed for periodic boundary conditions

  u0 =    u(j);    uM1 =    u(j-1);    uJp1 =    u(1);    uJp2 =    u(2)
  aVar0 = aVar(j); aVarM1 = aVar(j-1); aVarJp1 = aVar(1); aVarJp2 = aVar(2)

  !  this implicitly assumes that aVar is mostly positive, or else, mostly negative.
  if( aVar(1) > 0.0 ) then
     do i=j,3,-1
        u(i) = u(i) - (sigma/2) * ( 3*aVar(i)*u(i) - 4*aVar(i-1)*u(i-1) + aVar(i-2)*u(i-2) ) &
             + (sigma*sigma / 2 ) * ( 0.5*( aVar(i-1)+aVar(i) )*( aVar(i)*u(i)-aVar(i-1)*u(i-1) ) &
             - 0.5*( aVar(i-2)+aVar(i-1) )*( aVar(i-1)*u(i-1)-aVar(i-2)*u(i-2) )  )

        !  { aVar(i-1/2) * [aVar(i)u(i)-aVar(i-1)u(i-1)] - aVar(i-3/2) *[ aVar(i-1)u(i-1)-aVar(i-2)u(i-2) ]}
     enddo
     if ( iPeriodic .NE. 1) then
        u(2) = u(2) - sigma*( aVar(2)*u(2) - aVar(1)*u(1) )!  first order for u(2)
        u(1) = uLeft
     else ! periodic BC
        i = 2
        u(i) = u(i) - (sigma/2) * ( 3*aVar(i)*u(i) - 4*aVar(i-1)*u(i-1) + aVar0*u0 ) + (sigma*sigma / 2 ) * &
             ( 0.5*( aVar(i-1)+aVar(i) )*( aVar(i)*u(i)-aVar(i-1)*u(i-1) ) &
             - 0.5*( aVar0+aVar(i-1) )*( aVar(i-1)*u(i-1)-aVar0*u0 )  )
        i=1
        u(i) = u(i) - (sigma/2) * ( 3*aVar(i)*u(i) - 4*aVar0*u0 + aVarM1*uM1 ) + &
             (sigma*sigma / 2 ) * ( 0.5*( aVar0+aVar(i) )*( aVar(i)*u(i)-aVar0*u0 ) - &
             0.5*( aVarM1+aVar0 )*( aVar0*u0-aVarM1*uM1 )  )          
     endif
  else !  aVar(1) < 0 
     do i=1, j-2
        u(i) = u(i) + (sigma/2) * ( 3*aVar(i)*u(i) - 4*aVar(i+1)*u(i+1) + aVar(i+2)*u(i+2) ) &
             + (sigma*sigma / 2 ) * ( 0.5*( aVar(i+1)+aVar(i) )*( aVar(i)*u(i)-aVar(i+1)*u(i+1) ) &
             - 0.5*( aVar(i+2)+aVar(i+1) )*( aVar(i+1)*u(i+1)-aVar(i+2)*u(i+2) )  )
     enddo
     if( iPeriodic .NE. 1) then
        u(j-1) = u(j-1) - sigma*( aVar(j)*u(j) - aVar(j-1)*u(j-1) )!  first order for u(j-1)
        u(j) = uRight
     else ! periodic BC
        i = j - 1
        u(i) = u(i) + (sigma/2) * ( 3*aVar(i)*u(i) - 4*aVar(i+1)*u(i+1) + aVarJp1*uJp1 ) &
             + (sigma*sigma / 2 ) * ( 0.5*( aVar(i+1)+aVar(i) )*( aVar(i)*u(i)-aVar(i+1)*u(i+1) ) &
             - 0.5*( aVarJp1+aVar(i+1) )*( aVar(i+1)*u(i+1)-aVarJp1*uJp1 )  )
        i = j
        u(i) = u(i) + (sigma/2) * ( 3*aVar(i)*u(i) - 4*aVarJp1*uJp1 + aVarJp2*uJp2 ) &
             + (sigma*sigma / 2 ) * ( 0.5*( aVarJp1+aVar(i) )*( aVar(i)*u(i)-aVarJp1*uJp1 ) &
             - 0.5*( aVarJp2+aVarJp1 )*( aVarJp1*uJp1-aVarJp2*uJp2 )  )
     endif
  endif

  return
end subroutine UpWind2ndOrder_Variable_a


!-------------------------------------------------------------------------------------
! subroutine SubUpwind_NonCon_Form
!-------------------------------------------------------------------------------------


subroutine SubUpwind_NonCon_Form( u, dt, dz, aVar, j, uBottom,  uTop, iPeriodic, iUpWind, uTemp )

  ! as per communication from B. Wendroff on Jan 10, 2011, for:  du/dt =  - (d/dz)[ aVar(x) u ]
  ! see Word doc:   Advection in non uniform flow
  ! this version is similar to that in "linear lib" folder but source has been dropped

  implicit none

  integer j, i, iPeriodic, iUpWind
  real u(j), dt, dz, aVar(j), uBottom,  uTop, uTemp(j)
  real sigma
  real    u0,    uM1,    uJp1,    uJp2
  real aVar0, aVarM1, aVarJp1, aVarJp2

  ! Vol. 2 page 502 of Numerical Computation of Internal and External Flows by Hirsch
  ! du/dt = -aVar du/dz for aVAr constant
  ! sigma = aVar_max * dt / dz ,  sigma < 2 for stability if aVar is uniform

  ! u(i,n+1) = u(i) - (sigma/2)*[ 3*u(i) - 4*u(i-1) + u(i-2) ] + (sigma*sigma / 2 )*[ u(i) - 2*u(i-1) + u(i-2) ]

  !            all u's on RHS are at time step "n"

  if (j .le. 2) return

  sigma =     dt / dz

  u0 =    u(j);    uM1 =    u(j-1);    uJp1 =    u(1);    uJp2 =    u(2)  
  ! needed for periodic boundary conditions

  aVar0 = aVar(j); aVarM1 = aVar(j-1); aVarJp1 = aVar(1); aVarJp2 = aVar(2)  
  ! needed for periodic boundary conditions

  IF( iUpWind==1 ) THEN


     !  this implicitly assumes that aVar is mostly positive, or else, mostly negative.
     if( aVar(1) > 0.0 ) then
        do i=j,3,-1
           u(i) = u(i)      -   (sigma/2) *aVar(i) * ( 3*u(i) - 4*u(i-1) + u(i-2) )  &
                + (sigma*sigma / 2 ) * aVar(i)**2  * (   u(i) - 2*u(i-1) + u(i-2) )  &
                +     (sigma*sigma / 8 ) * aVar(i) * ( 3*u(i) - 4*u(i-1) + u(i-2) )  &
                * ( 3*aVar(i) - 4*aVar(i-1) + aVar(i-2) )             
        enddo
        if ( iPeriodic .NE. 1) then
           i = 2
           u(2) = u(2) - sigma * aVar(2) * ( u(2) - u(1) ) !  first order for u(2)
           u(1) = uBottom
        else ! periodic BC
           i = 2
           u(i) = u(i)    -   (sigma/2) *aVar(i) * ( 3*u(i) - 4*u(i-1) + u0 )  &
                + (sigma*sigma / 2 ) * aVar(i)**2  * (   u(i) - 2*u(i-1) + u0 )  &
                +     (sigma*sigma / 8 ) * aVar(i) * ( 3*u(i) - 4*u(i-1) + u0 )  &
                * ( 3*aVar(i) - 4*aVar(i-1) + aVar0 )  
           i=1
           u(i) = u(i)    -   (sigma/2) *aVar(i) * ( 3*u(i) - 4*u0 + uM1 )  &
                + (sigma*sigma / 2 ) * aVar(i)**2  * (  u(i)  - 2*u0 + uM1 )  &
                +     (sigma*sigma / 8 ) * aVar(i) * ( 3*u(i) - 4*u0 + uM1 )  &
                * ( 3*aVar(i) - 4*aVar0 + aVarM1 )   
        endif
     else !  aVar(1) < 0 
        do i=1, j-2
           u(i) = u(i) +   (sigma/2) *aVar(i) * ( 3*u(i) - 4*u(i+1) + u(i+2) )  &
                + (sigma*sigma / 2 ) * aVar(i)**2  * (   u(i) - 2*u(i+1) + u(i+2) )  &
                +     (sigma*sigma / 8 ) * aVar(i) * ( 3*u(i) - 4*u(i+1) + u(i+2) )  &
                * ( 3*aVar(i) - 4*aVar(i+1) + aVar(i+2) ) 
        enddo
        if( iPeriodic .NE. 1) then
           i = j - 1
           u(j-1) = u(j-1)-sigma * aVar(j-1)*( u(j) - u(j-1) )!  first order for u(j-1)
           u(j) = uTop
        else ! periodic BC
           i = j - 1
           u(i) = u(i)    +   (sigma/2) *aVar(i) * ( 3*u(i) - 4*u(i+1) + uJp1 )  &
                + (sigma*sigma / 2 ) * aVar(i)**2  * (   u(i) -2*u(i+1)  + uJp1 )  &
                +     (sigma*sigma / 8 ) * aVar(i) * ( 3*u(i) - 4*u(i+1) + uJp1 )  &
                * ( 3*aVar(i) - 4*aVar(i+1) + aVarJp1 )  
           i = j
           u(i) = u(i) +   (sigma/2) *aVar(i) * ( 3*u(i) - 4*uJp1 + uJp2 )  &
                + (sigma*sigma / 2 ) * aVar(i)**2  * ( u(i)   - 2*uJp1 + uJp2 )  &
                +     (sigma*sigma / 8 ) * aVar(i) * ( 3*u(i) - 4*uJp1 + uJp2 )  &
                * ( 3*aVar(i) - 4*aVarJp1 + aVarJp2 ) 
        endif
     endif

  ELSE  ! else use BW recommended central difference  ( modified by Source term )


     do i = 2, j-1
        uTemp(i) = u(i)  - (sigma/2) *aVar(i) * ( u(i+1) - u(i-1) ) &
             + (sigma*sigma / 4 ) * aVar(i) * &
             (  (aVar(i)+aVar(i+1))*(u(i+1)-u(i)) - (aVar(i)+aVar(i-1))*(u(i)-u(i-1))  )
     enddo
     if ( iPeriodic .NE. 1) then
        if( aVar(1) > 0 )then
           uTemp(1) = uBottom
           i=j
           uTemp(i) = u(i)  - sigma * aVar(i) * ( u(i) - u(i-1) )
        else
           uTemp(j) = uTop
           i = 1
           utemp(i) = u(i)  - sigma * aVar(i) * ( u(i+1) - u(i) )
        endif
     else ! periodic
        i=1
        uTemp(i) = u(i)  - (sigma/2) *aVar(i) * ( u(i+1) - u0 ) &
             + (sigma*sigma / 4 ) * aVar(i) * &
             (  (aVar(i)+aVar(i+1))*(u(i+1)-u(i)) - (aVar(i)+aVar0)*(u(i)-u0)  )
        i=j
        uTemp(i) = u(i)  - (sigma/2) *aVar(i) * ( uJp1 - u(i-1) ) &
             + (sigma*sigma / 4 ) * aVar(i) * &
             (  (aVar(i)+aVarJp1)*(uJp1-u(i)) - (aVar(i)+aVar(i-1))*(u(i)-u(i-1))  )
     endif   !  end of periodicity branch

     u(:) = uTemp(:)

  ENDIF  !  end of upwind versus centered choice

  return
end subroutine SubUpwind_NonCon_Form


!-------------------------------------------------------------------------------------
! subroutine NonLinearAdvection
!-------------------------------------------------------------------------------------

subroutine NonLinearAdvection( u, dt, dx, j, uLeft,  uRight, iPeriodic )

  ! see Word doc:   Advection in non uniform flow.doc
  !  /Users/harveyrose/Documents/MyLib/Advection in non uniform flow.doc
  !  advances one time step, dt,  du/dt = - u du/dx

  implicit none

  integer j, i, iPeriodic
  real u(j), dt, dx,  uLeft,  uRight
  real sigma
  real    u0,    uM1,    uJp1,    uJp2

  ! needed for periodic boundary conditions

  if (j .le. 2) return
  
  u0 =    u(j);    uM1 =    u(j-1);    uJp1 =    u(1);    uJp2 =    u(2)

  sigma =     dt / dx

  IF( u(1) > 0 ) THEN
     do i=j,3,-1
        u(i) = u(i) - 0.25E0 * sigma * ( 3 * u(i)**2 - 4 * u(i-1)**2 + u(i-2)**2 ) + &
             (1.E0 / 6.E0) * sigma * sigma * ( u(i-2)**3 + u(i)**3 - 2 * u(i-1)**3 )
     enddo
     if ( iPeriodic .NE. 1) then
        u(2) = u(2) - 0.5D0 * sigma * (  u(2)**2 - u(1)**2 )
        u(1) = uLeft
     else                        ! periodic and u > 0
        i=2
        u(i) = u(i) - 0.25E0 * sigma * ( 3 * u(i)**2 - 4 * u(i-1)**2 + u0**2 ) + &
             (1.E0 / 6.E0) * sigma * sigma * ( u0**3 + u(i)**3 - 2 * u(i-1)**3 )
        i=1
        u(i) = u(i) - 0.25E0 * sigma * ( 3 * u(i)**2 - 4 * u0**2 + uM1**2 ) + &
             (1.E0 / 6.E0) * sigma * sigma * ( uM1**3 + u(i)**3 - 2 * u0**3 )

     endif
  ELSE                             ! u < 0 
     do i=1, j-2
        u(i) = u(i) + 0.25E0 * sigma * ( 3 * u(i)**2 - 4 * u(i+1)**2 + u(i+2)**2 ) + &
             (1.E0 / 6.E0) * sigma * sigma * ( u(i+2)**3 + u(i)**3 - 2 * u(i+1)**3 )
     enddo
     if ( iPeriodic .NE. 1) then
        u(j-1) = u(j-1) - 0.5D0 * sigma * (  u(j)**2 - u(j-1)**2 )
        u(j) = uRight
     else                        ! periodic and u < 0
        i = j-1
        u(i) = u(i) + 0.25E0 * sigma * ( 3 * u(i)**2 - 4 * u(i+1)**2 + uJp1**2 ) + &
             (1.E0 / 6.E0) * sigma * sigma * ( uJp1**3 + u(i)**3 - 2 * u(i+1)**3 )
        i = j
        u(i) = u(i) + 0.25E0 * sigma * ( 3 * u(i)**2 - 4 * uJp1**2 + uJp2**2 ) + &
             (1.E0 / 6.E0) * sigma * sigma * ( uJp2**3 + u(i)**3 - 2 * uJp1**3 )
     endif
  ENDIF

  return
end subroutine NonLinearAdvection

!-------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------

